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PREFACE 


The  research  described  in  this  report  was  conducted  by  personnel 
of  the  Jacob  Blaustein  Institute  for  Desert  Research,  Ben-Gurion  University 
of  the  Negev,  Sede  Boqer,  Israel  from  1  July  1982  to  30  June  1983  under 
Grant  No.  AFOS  i8)285  to  the  Ben-Gurion  University  of  the  Negev,  Research 
and  Development  Authority,  P.0.  Box  1025,  Beersheva,  Israel. 

Participating  personnel  concerned  with  the  tasks  described  in  this  re¬ 
port  include  Prof.  Louis  Berkofsky,  Principal  Investigator,  Dr.  Avraham 
Zangvil,  Research  Associate,  Ms.  Ruth  Hildebrand,  Engineer-Programmer,  Ms. 
Andrea  Molod,  Meteorologist-Programmer,  and  Ms.  Perla  Druian,  Meteorologist- 
Programmer. 

Observational  data  used  in  this  study  were  obtained  from  the  Institute's 
micrometeorological  (4  m)  tower,  and  from  the  radiation  center,  collected  on 
a  data  logger  and  analyzed  in  the  laboratory.  Dust  data  were  obtained  from 
the  Institute's  Size  Selective  Inlet  High  Volume  Sampler. 

The  Director  of  the  Institute  during  the  conduct  of  this  study  was 
Prof.  Amos  Richmond. 


This  report  should  be  cited  as  follows: 

Berkofsky,  L.,  1983.  "The  Behayior  of  the  Atmosphere  In  the  Desert  Planetary 
Boundary  Layer,"  Final  Scientific  Report,  prepared  by  the  Desert  Meteorology 
Unit  of  the  Jacob  Blaustein  Institute  for  Desert  Research  for  the  U.S.  Air 
Force  Office  of  Scientific  Research,  AFSC,  Bolling  AFB,  D.C.,  20332. 
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Introduction 


In  the  previous  Final  Scientific  Report  (Berkofsky,  1982}  we  dis¬ 
cussed  the  physical  basis  for  the  formulation  of  the  Boundary  Layer  Model. 

In  the  present  report,  we  discuss  and  present  some  of  the  results  obtained 
from  various  versions  of  the  previous  model,  including  a  first  attempt  to 
include  dust  modelling.  We  then  re-derive  the  basic  equations,  using 
potential  temperature  instead  of  moist  static  energy  as  the  principal  thermo¬ 
dynamic  variable.  We  also  allow  the  lapse-rate  of  potential  temperature  to 
vary  with  time.  This  leads  to  the  necessity  for  a  prediction  equation  for 
the  lapse- rate,  which  is  derived. 

Some  Results  From  the  Initial  Model 
Numerical  Tests 

In  order  to  carry  out  numerical  integrations,  we  first  tested  a  simplified 
version  of  the  model  with  three  different  numerical  schemes.  The  equations 
we  used  were 

~  +  £>£  (1) 


-2- 


The  symbols  are  defined  in  the  previous  report.  These  are  very  similar 
to  the  "Shallow  Water"  equations.  We  tried  the  usual  "Leap-Frog"  scheme, 
which,  for  an  equation  of  the  form 


The  "Split-Explicit"  scheme  is  much  more  complicated  numerically,  but  essen¬ 
tially  consists  of  three  stages: 

a)  An  adjustment  stage,  in  which  only  the  time-dependent  geostrophic 
adjustment  terms  are  integrated  for  several  fractions  of  time  steps. 

b)  An  advective  stage,  in  which  only  the  time-dependent  advective  terms 
are  integrated  for  one  time-step,  using  the  results  of  a)  as  initial 
conditions. 
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c)  A  friction  stage,  in  which  only  the  time-dependent  friction  terms 
are  integrated  for  one  time-step,  using  the  results  of  b)  as  initial 
conditions. 

Our  results  show  that  there  is  very  little  difference  among  the  three 
schemes,  except  for  a  possible  slight  improvement  with  the  "split-explicit." 
However,  it  may  not  be  worth  the  extra  computer  time  to  continue  with  this 
method.  This  is  still  being  investigated.  If  we  abandon  the  "split-explicit," 
we  will  adopt  the  "Lax-Wendroff ,"  as  it  has  desirable  stability  characteristics. 


Dust  Concentration  Tests 

We  have  tested  a  simple  version  of  a  dust  concentration  model.  In  that 
version,  we  use  a  coupled  inversion  height  and  dust  concentration  equation, 
and  feed  in  whatever  other  variables  are  needed.  A  complete  description  of 
the  model,  and  the  derivation  of  the  inversion  height  and  dust  concentration 
equations,  is  given  in  the  attached  reprint  (Berkofsky,  1982). 

In  that  reprint,  we  derived  the  following  equation  for  the  variation  of 
inversion  height  with  time: 


ft.  C*.  foa  -~tV) 

U-fC)  +  (T6i-T4) 


(7) 


The  results  are  highly  dependent  upon  W^,  the  meso-scale  vertical  velocity, 
as  is  seen  from  Figs.  5  and  8  of  the  attached  paper.  There  are  several 
ways  of  resolving  the  problem  of  obtaining  suitable  values  of  Wh.  One  is 
to  interface  the  boundary  layer  model  with  a  free-air  mesoscale  model.  A 
second  method  is  based  on  the  following  considerations. 
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The  continuity  equation  when  used  in  Boussinesq  form,  i.e., 

V.  p  V  =  O  (8i 

may  lead  to  spurious  mass  losses  (Feliks  and  Huss,  1982).  To  avoid  this, 
the  continuity  equation  may  be  written 


V  p'V 


JL  Co- 


When  this  is  substituted  into  the  equation  for  ,  and  integrated  over 

the  depth  of  the  transition  layer  (h-k),  there  results 


^  -  04  \  +CurO  H  +  04  \ 


where 


(O0 


:^\  --  -1C  / 


*  A  formula  similar  to  this  was  derived  by  Lavoie  (1972),  but  in  that 
formula,  (W.  )u  is  absent. 


Thus,  in  this  form,  Wh  depends  upon  the  value  of  W  at  the  lower  boundary, 
(W^)g,  upon  the  convective  heating  and  divergence  of  flux  of  net  radiation 
(Wh)H,  and  upon  the  integral  of  the  horizontal  divergence,  (W^)^ .Preliminary 
estimates  indicate  that  (W^  is  the  dominant  term.  But  the  above  solution 
does  reveal  a  diurnal  variation  of  due  to  (W^,  the  heating  and  radiation 
parameters.  These  terms  may  conceivably  dominate  in  certain  situations.  We 
are  continuing  to  test  the  above  formulation  in  conjunction  with  Equation  (7). 

A  third  method  of  estimating  Wh  is  to  use  a  "single-station"  method 
(Tucker,  1973).  In  that  formulation,  and  using  our  modelling  assumptions 
from  the  previous  report,  we  find 


r«  +  Vr,  1 


All  terms  except  Q/  can  be  obtained  from  model  predictions.  However,  W.  is 

needed  to  calculate  .  This  may  be  done  interactively,  i.e.,  assume 

=  0.  Solve  for from  the  model  equation.  Put  this  into  the  Wh 

formula,  get  a  new  value,  re-sol ve  for  >  etc.,  until  convergence  is 

< it 

reached. 

An  alternative  to  the  above  is  simply  to  calculate  Wh  from  several  con- 


sectuve  radiosondes,  and  feed  this  information  into  the  model,  i.e., 


We  deal  with  the  system 


cm 


+  (uiaP)  -  ^  (Vur') 


It  +  +%(vlN)  +^fu'uj)  =-  ^('Vur') 


'S* 


(wc) 


The  symbols  are  as  previously  defined: 
u,  v,  w  =  components  of  the  velocity  vector. 
u',  V',  w'  =  turbulent  components  of  the  velocity  vector. 
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0  =  potential  temperature 
q  =  specific  humidity 
c  =  dust  concentration 

e%  q c'  =  turbulent  fluctuations  of  e,  q,  c.  Q(r)  is  the  velocity  of  c 
relative  to  the  air  at  any  one  height  and  time.  If  the  dust  concentration  at 
a  point  is  distributed  over  particles  of  different  sizes  and  fallspeeds, 
a(r)  at  the  point  is  then  a  kind  of  average. 


~  l  ^  F 


=  divergence  of  flux  of  net  radiation 


=  net  radiation  flux. 
The  bar  is  defined  by 


(*o 


=  S-  i -  °r  (21) 

where  is  any  scalar  variable.  The  horizontal  area  AXAy  is  large  enough 
to  contain  the  sub-grid  scale  phenomena,  but  small  enough  to  be  a  fraction 
of  the  mesoscale  system. 

We  integrate  the  above  system  (Eqs.  (14)-(19)  incl.)  through  the  in¬ 
version  layer,  of  thickness  £  above  the  inversion  layer  h.  Using  Leibniz's 
rule,  we  find  the  following  five  Interface  Conditions: 
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(  J  Au  4-  +4  =  -fu'ur) 


~  [AJ1^)  Air  4-  A  (u^  4-  4  (v~ =  -  fvW 


4-A(u&)^  ^  =-fur1 


*  M 


^  -HA  ( 


W  +  4 


=  -H 


— * 

“(ur-'f-xL^  AC  -v  -  4  (ut)  H  4  (v~c\  - 

\l  o>  '*v^ 

In  the  above  derivations,  we  have  assumed  that 

,4+f  . 

\  — X  -_LfF  -£)  = 

«U  or.  ^  T-^  v  U-t-i  W 


\  1>F 


(Cn.  4r  Ur')  „  .  r  =  P_a  -+  Ur') 


We  have  used  the  following: 


Parameterizations 

Winds 

Surface  Layer: 


U.fe')  =  ~-  A  (±+  h 


A 


,  0^-2^  (27) 


u*  =  friction  velocity,  kj  =  von  Kantian  constant 


Transition  Layer: 


Here 


!\  ( ■%')  la. 

Ifi')  V 


(28) 


k  ®  constant  height  of  surface  layer. 

Temperature 

(30) 

&(*,  f,%,o)zr  (+,  f,o\-hrfe)f?-4)  (31) 

=  f^-^+0-«-wt)]  r-k-4^  +.r(rts  02) 
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Moisture 


^'V*'0)  (%^.o)  +~s(0)( 2--k) 

*  ^*1  ■‘jt)  +  £?/»}  -  3  ft)]  (-t,--C)  +  ?m£ 

Dust 

C  h>  y  1  3 ,  t)  r  5  fe)  Cj,  /V,  ^ ,  t) 

=  [ptt+O-btt^c*  fM,P) 

Surface  Parameterizations 

(u’lv'  =  -  C*  4  ft)  1^1  ff  ) 
(v’ur'^-.Q  } 

(JVU  -  Cm.  A(4)u  (#„-  b) 

(WV  )*  =  Cte'^C' 


(urV  )<  -  Q  A-tt)  a  (c^-  Ck)j 


Cji  =  •  ooor  I  V\x  ,  \h 


In  the  above,  WgR  =  ground  soil  moisture,  eQR  =  ground  potential  temp¬ 
erature,  and  CGR  =  ground  dust  concentration,  must  be  predicted.  Wk  =  potential 
saturation  value  of  W. 


Ground  Albedo 


Rc-n  =  *  +  \> 


WK 


Ip  4o 


a,  b  constants 


Closures 


(uj'e1 )<  =-  a,(w-V \  (43] 

(.l/tf’C  =  ft*  if)  ^  (44] 

a . )  (45) 

We  now  apply  the  averaging  operator  Eq.  (29)  to  Equations  (14)-(19) 
inclusive,  introduce  the  five  Interface  Conditions  Equations  (22)-(26) 
Inclusive,  all  of  the  parameterizations  and  closure  assumptions,  Equations 
( 28) -( 30-45)  inclusive,  and  make  use  of  Leibniz's  rule,  and  obtain  the 
following  system  of  equations. 


uation 


Ground  Temperature  Equation 


r« 


=  -  >T‘’  H*  _  >2L  (g.  6  \ 

psc5a,  -  f.  l#i*  ®° 


where  represents  the  radiation  balance, 


Soil  Moisture  Equation 


(55) 


?  g- c.  (E±zll  _  cv(fw«-^) 


Tt 


(56) 


Pw  T, 

where  E^-P  represents  evaporation  minus  precipitation. 

In  deriving  the  above  system,  we  have  partitioned  the  radiative  fluxes 
and  the  moisture  fluxes  in  the  following  way:  the  actual  First  Law  of 
Thermodynamics  consisted  of  the  sum  of  Equations  (48)  and  (49).  We  then 
assumed  that  the  change  of  potential  temperature  at  the  top  of  the  surface 
layer  is  controlled  mainly  by  the  net  radiation  and  the  convective  heat 
flux  at  that  level,  while  the  change  of  lapse  rate  is  controlled  mainly 
by  the  net  radiation  and  the  convective  heat  flux  at  inversion  base  level. 

The  same  type  of  partition  led  to  Equations  (52)  and  (53). 

We  thus  have  10  prognostic  equations  in  the  10  unknowns  u,  v,  ek,  y  , 

(h-k) ,  qk,  c,  CR,  0gR,  WqR .  We  have  diagnostic  equations  for  calculating  fig, 

O  *  e 

g’  e* 


u 


A 


y\ 


(57) 
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We 


have  parameterized  I.U7  0  ,  (w  $  )o  ,  (,Ur  Q  ,  but  not  yet 


(W'V^-  We  can  obtain  F(k)  from  data,  but  we  have  not  specified 
yet  the  relation  between  F(h)  and  F(k),  if  any  (see  Kondratyev,  1965,  p. 

200).  RgR  is  parameterized.  We  also  need  Wk  and  CfiR.  p,  e2»  W2,  Wk 
are  prescribed.  A,  B,  A*,  D,  Cd,  CHQ,  ft  are  known.  The  lower  boundary 
conditions  have  been  discussed  in  the  previous  Final  Report. 

From  the  system  Equations  (46)-(56)  inclusive,  it  is  possible  to  set 
up  a  hierarchy  of  models,  as  described  in  the  previous  Final  Scientific 
Report.  Indeed,  the  results  of  the  dust  concentration  experiment  (Berkofsky, 
1982)  may  be  thought  of  as  a  zero  order  model. 

We  now  describe  a  second  version  in  the  hierarchy. 


A  Meso-Scale  Model 

A  considerable  simplification  of  the  basic  system  of  equations  results 
if  we  assume  no  wind  shear  through  the  inversion  layer.  Whether  or  not  this 
is  true  remains  to  be  verified.  This  assumption  implies  infinite  Richardson 
number,  making  the  inversion  layer  very  stable  with  respect  to  shear: 

We  also  assume  that  the  atmosphere  is  dry  and  that  there  is  no  wind 
shear  In  the  transition  layer.  There  is  actually  some  basis  for  the  latter 
assumption,  especially  during  daylight  hours  (Brown,  1974;  Martner  and 
Marwitz,  1982). 

With  these  assumptions,  the  equations  are 
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We  plan  to  Integrate  the  above  system  -  Equations  (57)-(65)  inclusive- 
over  a  horizontal  region  large  enough  to  encompass  a  substantial  fraction 
of  the  desert  in  this  region,  with  grid  spacing  as  yet  to  be  determined. 

The  lateral  boundary  conditions  will  be  of  the  "one-way  interactive"  type, 
i.e.,  the  data  from  any  outer  "coarse"  model  will  be  used  to  update  con¬ 
ditions  at  inflow  points  of  the  fine  model.  At  outflow  points  an  open 
boundary  condition  is  needed;  that  is,  a  boundary  condition  which  allows 
phenomena  generated  in  the  domain  of  Interest  to  pass  through  the  boundary 
without  undergoing  significant  distortion  and  without  influencing  the 
interior  solution.  Such  a  boundary  condition  is  the  Sommerfeld  radiation 
condition  (Orlanski,  1976;  Miller  and  Thorpe,  1981).  This  is  not  difficult 
to  apply  and  is  very  effective  In  removing  boundary  "noise"  and  reflection. 

We  shall  use  the  Lax-Wendroff  numerical  scheme,  given  by  Equations  (6) 
We  shall  apply  these  on  a  staggered  grid,  with  the  advective  form  of  the 
equations  at  the  half  time  steps,  and  the  flux  form  at  the  full  time  steps 
(Phillips,  1979). 

The  method  of  data  Initialization  has  not  yet  been  decided,  as  little 
has  been  developed  at  this  point  for  desert  and  tropical  regions. 
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Conclusions 


We  have  re-derived  a  general  system  of  vertically  integrated  equations 
for  the  desert  planetary  boundary  layer,  including  an  equation  for  prediction 
of  dust  concentration.  Innovative  features  of  the  revised  model  are  the 
derivation  of  equations  for  prediction  of  lapse  rate  of  temperature  and  moisture, 
respectively. 

A  test  of  a  zero  order  model,  in  which  inversion  height  and  dust  concen¬ 
tration  are  predicted,  while  other  data  are  fed  in,  leads  to  what  appear  to 
be  realistic  results.  These  results  were  highly  dependent  upon  the  time  vari¬ 
ation  of  meso-scale  vertical  velocity,  which  was  prescribed.  Several  methods 
for  obtaining  this  vertical  velocity  from  observed  data  are  proposed. 

We  have  derived  another  system  in  the  model  hierarchy.  This  system  is 
two-dimensional  in  space,  is  representative  of  a  dry  atmosphere,  with  no  wind 
shear  across  the  inversion  layer.  This  system  will  be  integrated  over  a 
limited  area  of  the  desert,  using  a  fine  mesh  with  Sommerfeld  radiation  out¬ 
flow  conditions.  Data  initialization  techniques  peculiar  to  the  desert  will 
be  used. 

Recommendations 

We  have  gathered  some  ground-based  dust  data,  using  a  Sierra  Instruments 
Size  Selective  Inlet  High  Volume  Cascade  Impactor.  After  much  difficulty 
with  shredding  of  the  filter  paper  (substrates),  we  discovered  that  teflon 
backed  substrates  are  now  produced  in  the  U.S.,  and  we  plan  to  continue  our 
dust  measurements  using  these. 


We  are  taking  radiosonde  and  pilot  balloon  observations,  and  plan 
to  put  into  operation  a  tethersonde.  We  operate  radiation  devices  on  a 
regular  basis,  and  do  intermittent  heat  budget  studies.  All  of  the  above 
will  be  used  to  determine  various  parameters  and  profiles  for  the  model 
being  tested. 

We  expect  to  interact  with  other  people  in  Israel  to  study  the  effects 
of  meteorology  on  transmission  of  radiation  and  imaging  through  natural 
dusty  environments. 
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ABSTRACT 

A  boundary  layer.model,  to  investigate  various  aspects  of  desert  circulation,  has  been  devised.  The  planetary 
boundary  layer  is  divided  into  a  surface  (constant  flux)  layer,  a  transition  layer  and  an  inversion  layer.  The 
basic  equations  (motion,  continuity,  first  law  of  thermodynamics,  dust  concentration  equation)  are  integrated 
with  respect  to  height  from  the  bottom  of  the  transition  layer  up  to  the  top  of  the  boundary  layer.  It  is 
assumed  that  this  height  coincides  with  that  of  the  inversion,  which  is  allowed  to  vary  in  space  and  time. 
Interface  conditions  are  derived  by  integration  of  the  equations  across  the  inversion  layer.  These  may  be 
used  to  eliminate  the  momentum,  moisture  and  dust  concentration  turbulent  fluxes  from  the  vertically 
integrated  equations.  The  turbulent  heat  flux  is  eliminated  from  the  first  law  by  a  closure  assumption,  relating 
the  turbulent  flux  at  the  base  of  the  inversion  to  that  at  the  bottom  of  the  transition  layer.  Modelling 
assumptions,  specifying  the  variation  with  height  of  wind,  temperature  and  dust  concentration,  are  used  in 
the  vertically  integrated  equations.  The  resulting  equations  predict  the  horizontal  wind  components  at  a 
“mean”  level,  and  the  temperature  and  dust  concentration  at  the  bottom  of  the  transition  layer. 

A  simplified  version  of  the  model,  in  which  advective  terms  are  neglected,  couples  the  dust  concentration 
and  inversion  height  equations.  In  the  former,  the  time  variation  of  the  dust  concentration  depends  upon 
the  difference  between  sedimentation  and  turbulent  diffusion  of  dust,  and  upon  the  time  variation  of  the 
inversion  height.  The  latter  depends  upon  the  turbulent  transport  of  heat  from  the  ground,  upon  the  difference 
between  the  initial  temperature  and  temperature  at  a  later  time,  and  upon  the  vertical  velocity  of  the  larger 
scale  motions  at  interface  height  A  numerical  integration  of  this  system  describes  the  evolution  of  the 
inversion  height  its  effect  on  the  dust  concentration  evolution,  at  well  as  the  evolution  of  the  dust  concen¬ 
tration  itself. 


Introduction 

The  problem  of  dust  transport  has  been  receiving 
icreasing  attention  in  recent  years.  Of  particular  in¬ 
rest  has  been  the  transport  of  Saharan  dust  by  large- 
ale  atmospheric  circulations  (Morales,  1979). 
Smaller  scale  transport  is  also  of  interest  for  a  va¬ 
sty  of  reasons.  Among  these  are  the  siting  of  build- 
gs,  location  of  agricultural  areas,  effects  on  health, 
Fects  on  aviation,  effects  on  machinery,  and  effects 
i  light  transmission  by  various  atmospheric  probes. 
To  model  several  scales  simultaneously  is  a  for- 
iidabie  task.  The  few  studies  that  have  been  made 
-rnfine  themselves  to  a  single  scale.  Since  we  are  in- 
rested  in  the  dust  problem  as  it  concerns  future 
ttiement  of  the  Negev  (Berkofsky,  et  al.,  1981),  we 
’  all  confine  ourselves  to  the  mesoscale,  commen- 
irate  with  the  size  of  the  Negev. 

The  Negev  is  under  the  influence  of  the  descending 
■ancfa  of  a  subtropical  high  pressure  cell.  As  a  result, 
tbsidence  inversions  exist  on  the  average  222  days 
sr  year  (Shaia  and  Jaflfe,  1976).  Even  when  no  dust 
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storm  influences  the  region,  the  loess  which  is  picked 
up  by  the  wind  is  trapped  within  the  boundary  layer. 
To  be  able  to  predict  the  daily  variation  of  dust  con¬ 
centration  requires  a  method  to  predict  the  variation 
of  inversion  height. 


2.  The  model 

Although  we  have  developed  a  model  which  is 
completely  time  and  space  dependent,  we  shall  deal 
here  with  a  very  simple  version,  mainly  as  a  testing 
ground  for  the  complete  model. 

The  method  consists  of  averaging  vertically  the 
basic  equations,  and  introducing  modelling  assump¬ 
tions  in  the  vertical,  thereby  eliminating  the  vertical 
dimension.  At  the  same  time,  various  parameteriza- 
tions  are  introduced  for  the  fluxes  of  dust  and  heat. 

We  divide  the  atmospheric  boundary  layer  into 
layers,  as  shown  in  Fig.  1.  Within  the  surface  layer 
of  height  z  -  k,  we  assume  that  the  wind  obeys  the 
constant  flux  formula.  The  bright  of  the  transition 
layer  is  h  -  k.  The  thickness  of  the  inversion  layer 
is  6. 

The  temperature  variations  above  the  surface  layer 


Flo.  1.  Schematic  of  planetary  boundary  layer. 


are  defined  by 

T\z,  t)  -  Tk(t)  +  -tfzXz  “  k),  (l) 

T\z,  0)  -  7*/  +  y,(z  -  k),  (2) 

where  7,,  a  constant,  is  the  increase  of  temperature 
with  height  at  t  *  0  and  also  above  base  of  inversion 
at  all  times,  y(z)  is  the  falloff  of  temperature  with 
height  in  transition  layer,  Tk{t)  the  temperature  at 
z  -  k,  and  Tu  the  temperature  at  height  z  -  k  at 
initial  time. 

We  then  find  that  the  inversion  strength  is  given 
by 

AT  «  ( Tu  -  T*)  +  [7.  -  t(A)P  -  k)  +  t (3) 
In  terms  of  potential  temperature  6,  we  have 

AB  -  <«w  -  **)  +  (^)'l7i  -  7</i)P  -  *) 

+  (4) 

where  p  is  pressure,  x  »  /I  is  the  gas  constant, 
and  e,  the  specific  heat  at  constant  pressure. 

The  variation  of  dust  concentration  with  height  is 
given  by 

Qz,  t )  -  D(z)Ck{t),  (5) 

D(z)  -  (z/kY«*>V'\  (6) 

Hoe  kt  is  the  von  Kirmin  constant  (0.4),  £/,  the 
friction  velocity,  and  Q(r)  the  sedimentation  or  set¬ 
tling  velocity  of  dust  particles  of  radius  r. 

Gillette  and  Goodwin  (1974)  strived  a  one-dimen¬ 
sional  steady-state  equation  in  which  vertical  turbu¬ 
lent  diffusion  of  dust  balanced  sedimentation  in  the 
surface  layer.  Their  solution  was  of  the  form  shown 
as  Eq.  (5),  with  the  time-dependence  omitted,  and  in 


which  D(z)  was  given  by  (6),  except  that  their  result 
used  Zo,  the  roughness  height,  instead  of  k.  We  have 
assumed  that  a  time-dependent  relation  of  this  form 
exists  throughout  the  transition  layer.  The  form  of 
D(z)  should  be  determined  from  data,  which  are  at 
present  unavailable.  For  the  time  being  we  use  the 
form  above,  as  a  modelling  assumption. 

The  equation  for  the  rate  of  change  of  dust  con¬ 
centration  in  a  vertical  column  is 

£+£ic,*+«,i--£(^  (7) 

where  w  is  the  vertical  velocity  of  the  air.  The  velocity 
of  C  (the  dust  content  oi  the  air)  relative  to  the  air 
at  any  one  height  and  time  is  represented  by  0 (r),  a 
negative  quantity  since  it  is  downward.  If  the  dust 
concentration  at  a  point  is  distributed  over  particles 
of  different  sizes  and  fallspeeds,  fl(r )  at  that  point  is 
then  a  kind  of  average,  w'  and  C'  are  deviations  from 
their  means,  and  the  bar  is  a  horizontal  average  over 
a  sufficiently  huge  area. 

Eq.  (7),  it  is  noted,  retains  the  vertical  motion 
term,  but  omits  the  horizontal  advection.  Thus,  the 
omisaon  of  V  •  (CV)  -  CV  •  V  +  V  •  VC  and  V  •  (gV) 
*  gV*  V  +  V*  V#  in  the  first  law  of  thermodynamics, 
to  be  considered  later  (here  V  is  the  horizontal  wind 
vector),  is  tantamount  to  assuming  that  advection 
and  horizontal  divergence  are  in  balance.  The  hori¬ 
zontal  divergence  need  not  be  zero,  so  that  the  vertical 
motion  term  need  not  be  neglected.  Our  aim,  at  this 
stage,  is  not  to  develop  an  operational  model,  but  to 
test  the  feasibility  of  doing  so.  In  steady-state  models, 
like  those  of  Gillette  and  Goodwin  (1974)  the  vertical 
motion  term  is  correctly  omitted,  since  they  assumed 
zero  horizontal  divergence.  Here  we  are  interested  in 
seeing  the  effect  of  vertical  motion  on  the  evolution 
of  inversion  height  and  dust  concentration,  without 
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having  to  solve  the  complete  set  of  primitive  equa¬ 
tions.  For  this  reason,  we  include  the  vertical  motion 
term  in  the  dust  concentration  and  in  the  first  law  of 
thermodynamics.  For  purposes  of  computation  in  the 
final  equations,  we  shall  assume  a  distribution  with 
time  of  vertical  velocity  at  inversion  height 


We  apply  the  operator 

*  ,8> 

to  Eq.  (7),  and  introduce  Eq.  (5).  Integrating  by  parts 
using  Leibniz's  rule,  and  recalling  that  h  -  hit),  yields 


dC* 

dt 


CAD(h)  -D]  j(  +  [(w*  +  0)  -  D(h\w„  +  Q )]C*  -  (w'C%  +  (w'O* 

_____ 


(9) 


We  parameterize  (*v'C)k  by  means  of 

-(wV)  -  A*(z)  — ,  (10) 

i.c.> 

~WC%  »  A  *(hj^)Ck  .  00 

where  A*  (z)  is  the  coefficient  of  eddy  diffusion,  given 
by 

A*(z)  «  Umkyz ,  z  <  k 

[see  Gillette  and  Goodwin  (1974)]  '  (12) 

A*(z)  =  constant  k<z<h 

and  (w'CO*  by  means  of 


(w'C%  -  Cw«*(Cor  -  C*).  (13) 

Hoe  Con  is  the  concentration  in  a  thin  layer  near  the 
ground  and  uk  is  the  east-west  component  of  wind 
at  z  -  k.  This  parameterization  for  dust  concentra¬ 
tion,  analogous  to  that  for  momentum  and  heat,  pro¬ 
vides  the  mechanism  for  vertical  turbulent  transport 
of  dust  In  (13),  Crfis  the  transfer  coefficient  for  dust. 
For  sufficiently  small  particles  it  may  be  assumed  that 
Gris  the  same  as  the  transfer  coefficient  for  momen¬ 
tum  (Cbepil  and  Woodruff,  1957).  Then 

Gr - 0.0005  Vm,  Kin  ms"1.  (14) 

Then,  finally 


sedimentation 


diffusion 


[D(h)  -  D)Ck  —  +  {[(w*  +  fl)  -  £>(AX"a  +  0))  +  A*D\h))Ck  +  QMGm  -  Ck) 

D(h-k)  ' 


The  terms  responsible  for  sedimentation  and  diffu¬ 
sion  are  labelled.  Notice  that  the  time  rate  of  change 
of  inversion  height  has  an  effect  on  the  concentration 
tendency.  But  this  effect  is  small.  However,  the  quan¬ 
tity  h  -  k  in  the  denominator  strongly  influences 
dCJdt  since  h  -  k  goes  through  a  large  variation  dur¬ 
ing  the  course  of  a  day. 

The  sedimentation  velocity  0(r)  depends  upon  the 
radius  of  the  dust  particle.  Since  the  above  equation 
refers  to  a  single  value  of  Q,  it  predicts  the  concen¬ 
tration  as  if  there  were  a  single  particle  size  within 
the  distribution.  To  account  for  other  sizes  would 
require  solving  the  equation  for  other  values  of  Q(r ), 
and  summing  the  results.  Of  course,  initial  values  of 
Ck  would  have  to  be  known  for  various  size  particles. 
If  we  are  concerned  with  the  case  of  atmospheric  dust 
particles  (see  Fig.  2  (after  Bagnold,  1954)],  we  see 
that  these  are  sufficiently  small  so  that  they  have  only 
a  slight  sedimentation  velocity.  Thus,  what  is  likely 
is  that  the  heavier  particles  ffill  out,  while  the  smaller 
particles,  which  have  smaller  sedimentation  veloci¬ 
ties,  require  only  a  slight  upward  motion  or  upward 
transfer  to  keep  them  in  suspension. 


The  derivation  of  an  appropriate  inversion  height 
equation  has  been  the  subject  of  intensive  investi¬ 
gation  for  many  yean  (see,  for  example,  Tennekes 
and  Driedonks,  1981;Smeda,  1979;  Stull,  1973).  Ten¬ 
nekes  (1973),  Betts  (1973)  and  Carson  (1973)  have 
assumed  the  entrainment  rate  of  heat  from  the  in¬ 
version  layer  is  such  that 

(16) 

where  (wW),  in  the  surface  value.  We  shall  use  this 
relation  in  the  form 

-Wk-A&Tk.  (17) 

We  make  use  of  this  relation  in  the  “interface  con¬ 
dition”  derived  in  the  following  way.  The  first  law  of 
thermodynamics  can  be  written 

£+£<*»- ns) 

where  Q  is  the  divergence  of  the  flux  of  net  radiation. 
We  integrate  this  equation  between  the  levels  h  and 
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Bo.  2.  Relative  size  and  fall  rate  (m  a'1,  ordinate  of  lower  figure)  of 
small  particles  (after  Bagnniri,  1954). 


h  +  i,  and  let  5  —  0.  Then 


(ft  -  w*)a0  «  -(w'0%  •  (19) 

Using 

0 W)k  -  CHoUttfoa  -  «*),  (20) 

where  foa  is  the  potential  temperature  at  the  ground, 
and  Cho  the  heat  transfer  coefficient,  together  with 
Eqs.  (19)  and  (4),  we  obtain 

dh 

—  m  w* 

at 

•4|Cho«*(®ob  _  8*) 

te)lY' 

-  yihWh  -  k)  +  (9U  -  0k)  + 

(21) 

which  is  similar  to  the  Tennekes  (1973)  model. 

Eqs.  (21)  and  (15)  are  two  of  the  equations  of  the 
complete  prediction  scheme.  To  compile  the  system, 
we  would  need  five  more  equations:  the  first  equation 
of  motion  for  prediction  of  w*,  the  first  law  of  ther¬ 
modynamics  for  prediction  of  7*,  a  ground  temper¬ 
ature  equation  for  prediction  of  7o*,  the  equation  of 
continuity  for  calculation  of  w*,  and  an  equation  for 
prediction  of  dust  concentration  at  the  ground  Can. 

We  have,  in  fact,  derived  this  system,  using  the 
same  approach  as  above.  Here  we  wish  to  discover 
whether  it  is  possible  to  obtain  meaningful  results 
from  an  abbreviated  version,  before  going  ahead  with 
the  frill  model  Thus,  wherever  possible,  we  shall  feed 
in  observed  data,  such  as  u*.  Tan .  7V  The  problem 


of  obtaining  observed  data  for  COR  is  more  compli¬ 
cated,  and  more  property  belongs  to  the  realm  of  sed- 
imentology.  However,  in  discussions  with  Prof.  P.  R. 
Owen  (1981,  private  communication),  the  assump¬ 
tion  of  setting  CaR  -  constant  for  24  h  is  reasonable. 
There  is  obviously  much  room  for  improvement  of 
this  assumption,  and  we  do  not  make  any  strong 
claims  for  it  at  this  point. 


Table  I.  Oita  taken  at  Sede  Boqcr,  Israel,  18  June  1979. 
The  Wk  values  are  assumed. 


time 

(LSD 

(7w  -  7*) 

ro 

MTok  -  7*) 
(cm  s*’  *C) 

(cm  s'1) 

w* 

(cm  s'1) 

0900 

0.00 

850 

291 

1.0 

1000 

-1.95 

4738 

480 

1.0 

1100 

-3.35 

5440 

515 

1.5 

1200 

-4.35 

5566 

550 

1.75 

1300 

-5.75 

5470 

560 

2.0 

1400 

-5.05 

5220 

570 

2.0 

1500 

-4.75 

3540 

695 

2.0 

1600 

-4.05 

2360 

820 

1.0 

1700 

-2.50 

1320 

680 

0.0 

1800 

-.89 

760 

540 

-1.0 

1900 

.95 

340 

418 

-1.0 

2000 

2.65 

20 

295 

-1.0 

2100 

4.00 

-220 

223 

-1.5 

2200 

5.27 

-330 

151 

-2.0 

2300 

6.35 

-80 

204 

-2.0 

0000 

7.25 

18 

258 

-2.0 

0100 

8.05 

-200 

223.5 

-2.0 

0200 

8.88 

-580 

189 

-2.0 

0300 

9.95 

-590 

190 

-1.0 

0400 

10.55 

-535 

192 

-1.0 

0500 

10.45 

-400 

166 

-1.0 

0600 

9.58 

-220 

140 

0.0 

0700 

6.55 

80 

120 

0.25 

0100 

2.70 

570 

101 

0.50 

r 

i 
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Since  (Po/PkY  <*»  1  for  a  range  of  values  of  h,  we 
use  that  here.  Further,  since  (Pq/P*Y  *»  1  and  ytS 
<  (fiu  -  0*)  for  reasonable  values  of  both  these  quan¬ 
tities,  we  may  write  (21)  as 

dh  A,CHoUt(Tax  -  7a) 

[yx-yVMh-k)  +  (Tu-Tk)'  K  1 

where  Ton  is  ground  temperature. 

Eqs.  (22)  and  (15)  comprise  the  test  scheme,  given 
certain  input  variables  and  constants. 

3.  Calculations  and  results 

Initial  values  of  A*  300 m  at  0900  LST  and 
Ck  “  100  MS  no’3  are  assumed.  The  latter  value  is 
taken  from  Junge(  1979)  as  representative  of  Saharan 
air.  We  use  A,  *  0.2,  y,  **  0.5  X  10-4  K  cm"',  y(h) 
-  -0.5  X  10’4  K  cnT\  C«o  -  1-695  X  10’\  A*fh) 
■  1.6  x  104  cm2  s’1,  Cor  »  100  pg  m_J  -  constant, 
k  *  20  m.  We  also  assume  wk  -  0.1  w*.  Table  1  gives 
a  tabulation  of  Tk,  -  Tk,  uk(TQ%  -  Tk)\  uk,  and  an 
assumed  distribution  h»*  (Fig.  3)  for  Sede  Boqer,  Israel 
for  1 8  June  1 979.  The  value  of  A*(h)  is  obtained  from 
Eq.  (12)  with  £7*  -  20  cm  s"\  kx  *  0.4,  r  *  20  m 
(top  of  surface  layer). 

The  distribution  of  wk  is  problematical.  These  val¬ 
ues  should  come  (as  stated  earlier)  from  solutions  of 
the  continuity  equation.  For  sensitivity  test  purposes, 
we  have  assumed  rising  motion  during  daylight  hours, 
and  sinking  motion  from  late  afternoon  until  sunrise. 
This  is  certainly  qualitatively  correct  on  a  clear  sum¬ 
mer  day,  but  foe  values  could  be  changed  somewhat 
without  greatly  affecting  the  shapes  of  foe  calculated 
curves  for  h  and  C*.  The  notable  difference  occurs 
when  w*  -  0  at  all  times. 

We  have  carried  out  calculations  for  three  different 
values  of  0(r):  -1.0,  -0.1  and  -0.01  cm  s"'.  As  can 
be  seen  from  Fig.  2,  these  represent  particles  of  radii 
10"*,  10’3  and  ♦  x  10-4  mm,  i,e.,  approximately  10, 
1  and  0.4  pm — all  within  foe  range  of  atmospheric 
dust  particles.  Fig.  4  (after  Gillette)  shows  a  graph  of 
QJU,  vs.  foe  ratio  of  upward  fluctuations  to  down¬ 
ward  fluctuations.  C/«  is  a  measure  of  the  standard 
deviation  of  w'.  We  see,  for  foe  values  used  here 


(Um  =  20  cm  s'*.  0  =  -1.0,  -0.1,  -0.01  cm  s'1),  that 
the  ratio  of  upward  to  downward  fluctuations  is  very 
close  to  1.  This  means  that  particles  with  these  radii 
will  probably  stay  in  suspension  unless  acted  on  by 
sufficiently  large  vertical  velocity  of  the  larger  scale 
motions. 

Fig.  5  shows  foe  results  of  solving  Eq.  (2 1 )  with  the 
data  given.  On  foe  same  figure  is  the  march  of 
(p*c,r*  X  Heat  Flux  (this  curve  is  labelled  HF), 
near  foe  ground.  Although  no  verification  data  are 
as  yet  available,  it  can  be  seen  that  the  variation  of 
inversion  height  throughout  the  day  is  reasonable. 
The  maximum  of  1051  m  is  reached  at  1800  LST. 
This  maximum  lags  that  of  the  heat  flux  by  ~6  h, 
which  is  in  agreement  with  the  O’Neill,  Nebraska  re¬ 
sults  shown  in  Fig.  1  of  Canon  (1973).  Our  results 
are  very  sensitive  to  vertical  velocity  variations  (see 
Figs.  6-9). 

Fig.  6  shows  the  results  of  solving  Eq.  (22)  together 
with  (15),  for  fl(r)  *  -1.0  cm  s’1  and  0(f)  =  -0-1 
cm  s-1.  The  concentration  Ck  decreases  as  the  day 


Flo.  4.  Sedimentation  velocity  compared  to  vertical  velocity 
fluctuation*  (after  Gillette,  1979). 
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Tim*  (hr) 

Fig.  5.  Time  variation  of  inversion  height  (ft )  and  surface  heat  flux  (HF). 

progresses,  reaches  a  minimum  around  2  h  before  the  time,  there  is  a  weak  upward  flux  of  C,  i.e.,  Cgr 
maximum  of  inversion  height,  and  then  goes  back  to  >  C*,  to  increase  Ck.  Also  h  is  decreasing,  so  that  the 
a  slightly  higher  value  than  the  initial  one  by  the  end  dust  is  being  “squeezed”  into  a  smaller  volume.  No- 
of  the  24  h  period.  What  appears  to  be  happening  is  tice  that,  for  Q  -  -0.1  cm  s"\  Ck  stays  higher  than 
that  when  (w*  +  0)  >  0  (from  0900-1600),  the  par-  for  0  -  -1.0,  simply  because  of  the  smaller  settling 
tides  rise.  C*  decreases  during  that  time  because  the  velocity.  The  results  for  0  =  -0.01  are  almost  the 
inversion  height  increases.  Even  though  there  is  an  same  as  for  Q  -  -0.1. 

upward  transfer  of  dust  (C**  >  Ct),  the  increase  of  Fig.  7  shows  the  total  concentration  f}  Cdz  for  the 
inversion  height  phis  the  net  upward  motion  distrib-  three  values  of  Q.  In  all  cases,  the  total  concentration 
utes  the  dust  upward.  After  1600,  when  (w*  +  Q)  increases  and  then  decreases.  This  is  undoubtedly  due 
<  0  the  particles  sink,  and  Ck  increases.  At  the  same  to  the  fact  that  we  are  forcing  COR  to  remain  constant. 


Tim*  (hr) 


Fta.  6.  Tim*  variation  of  Invenion  height  (ft)  and  dust  concentration  (Ct) 
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Tim*  (hr) 

Flo.  7.  Time  variation  of  total  concentration. 


In  actual  fact,  as  long  as  there  is  an  upward  transfer  The  lighter  particles  are  distributed  upward  more 
(Cn  >  Ck),  Cqr  should  decrease.  The  fact  that  we  easily,  increase  the  total  concentration,  and  stay  in 
are  keeping  CGR  constant  is  tantamount  to  assuming  suspension  longer.  This  is  so  because  the  subsidence 
that  a  source  exists  to  maintain  the  value  of  Cor-  velocity  plus  the  sedimentation  velocity  of  these  par- 
Thus  dust  is  being  added,  and  the  total  concentration  tides  is  less  than  for  the  heavier  particles, 
increases.  Later  in  the  day,  the  subsidence  overcomes  Fig.  8  shows  the  results  when  w*  *  0  at  all  times, 
the  upward  turbulent  transport,  and  decreases  the  The  inversion  height  h  reaches  a  somewhat  smaller 
total  concentration.  To  account  for  the  changes  in  maximum  than  when  wh  +  0,  at  about  the  same  time. 
Cor  requires  a  more  accurate  prediction  equation  and  remains  constant  after  that.  The  concentration 
than  simply  assuming  CGr  to  be  a  constant  C*  decreases  steadily  with  time  when  ft  =  ~  1.0  cm 


Tims  (hr) 

variation  of  invcrtion  height  (A)  and  dust  concentration  (C*)(m>«  *  0). 
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s'1,  but  remains  constant  when  0  =  -0.1-0.01  cm 
s'1.  This  simply  means  that  the  lighter  particles  are 
essentially  in  suspension  when  not  acted  on  by  sig¬ 
nificant  vertical  motion.  The  role  of  subsidence  in 
determining  inversion  height  and  concentration  is 
clearly  shown  by  these  results. 

Fig.  9  shows  the  total  concentration  for  the  above 
three  values  of  Q  when  w*  «  0.  Similar  arguments  to 
the  above  can  be  made  to  explain  the  results. 

To  date,  no  data  exist  at  Sede  Boqer  to  verify  any 
of  the  above  results. 


4.  Coactaskms 

An  extremely  simplified  version  of  a  model  for 
prediction  of  dust  concentration  has  yielded  what 
appear  to  be  reasonable  results.  We  are  thus  encour¬ 
aged  to  continue  with  the  full  model,  in  which  the 
entire  interactive  system  (primitive  equations)  will  be 
used.  Consideration  will  also  be  given  to  possible 
heating  effects  of  lower  troposphere  dust  concentra¬ 
tions. 

At  the  same  time,  we  are  developing  a  program  for 
dust  and  inversion  height  measurements  at  our  lo¬ 
cation  in  the  Negev.  We  will  then  be  able  to  verify 
present  and  future  results,  and  use  the  data  for  other 
dust  studies,  as  mentioned  in  the  Introduction. 

Acknowledgment.  1  wish  to  thank  Mr.  Eli  Zangvil 
far  programming  the  model  and  carrying  out  the  cal¬ 
culations. 
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